MOwNiT - laboratorium 8

Całkowanie

https://github.com/kzajac/Mownit2018/blob/master/Mownit_lab8.ipynb

In [1]:
using QuadGK
using Polynomials
using Plots
using Roots

Wykład

Kwadratury Gaussa bazują na tzw. odciętych punktów Gaussa $x_{i}$ oraz współczynnikach $a_{i}$

In [2]:
# funkcja obliczająca odcięte punktów Gaussa oraz ich współczynniki 
(xp,a)=gauss(Float64,7)
Out[2]:
([-0.949108, -0.741531, -0.405845, 0.0, 0.405845, 0.741531, 0.949108], [0.129485, 0.279705, 0.38183, 0.417959, 0.38183, 0.279705, 0.129485])

Podstawowe użycie w/w danych: aby policzyć $\int_{-1}^{1} f dx$ używamy odciętych punktów wraz z wagami wg wzoru $\sum_{i=1}^{n}{a_{i}f(x_{i})}$

In [3]:
f(x)=x^2
sum(a .* f.(xp))
Out[3]:
0.6666666666666669

Zadanie 1

  • Korzytając z pakietu Polynomials zaimplentuj wielomiany Legendre'a zdefiniowane w taki sposób:

$P_{0}(x)=1$

$P_{1}(x)=x$

$P_{k+1}(x)=\frac{2k+1}{k+1} x P_{k}(x)-\frac{k}{k+1}P_{k-1}(x)$

  • narysuj ich wykresy w przedziale (-1,1).
  • sprawdź(np. za pomocą fukcji roots z pakietu Polynomials), że ich zera sa odciętymi punktów Gaussa. Test wykonaj dla wielomianow od 2 do 4 stopnia
  • podaj związek tego faktu z podstawowym twierdzeniem kwadratur Gaussa (z wykładu)
In [4]:
function legendre(n)
    p0 = Poly([1.])
    p1 = Poly([0., 1.])  
    if n == 0
        return p0
    end
    for k = 1:(n - 1)
        p = Poly([0., (2 * k + 1.)]) * p1
        p2 = (p - (k * p0)) / (k + 1.)
        p0 = p1
        p1 = p2
    end
    return p1
end
Out[4]:
legendre (generic function with 1 method)
In [5]:
x = -1.:0.05:1.
for i = 2:5
    y = legendre(i)(x)
    p = plot!(x, y, label = "P_$i")
end
plot!(size = (950, 500), legend = :bottomright)
Out[5]:
-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 P_2 P_3 P_4 P_5
In [6]:
for i = 2:4
    println("degree: ", i)
    r = find_zeros(legendre(i), -1., 1.)
    println("\troots: ", r)
    (xp, a) = gauss(Float64, i)
    println("\tgauss: ", xp)
end

degree: 2
	roots: [-0.57735, 0.57735]
	gauss: [-0.57735, 0.57735]
degree: 3
	roots: [-0.774597, 0.0, 0.774597]
	gauss: [-0.774597, 0.0, 0.774597]
degree: 4
	roots: [-0.861136, -0.339981, 0.339981, 0.861136]
	gauss: [-0.861136, -0.339981, 0.339981, 0.861136]

Zadanie 2

  • Napisz funkcję, która dla podanej liczby punktów Gaussa $k$ oraz funkcji $f$ policzy $\int_{-1}^{1} f dx$ metodą używającą funkcji:
    gauss(k)
    oraz sumy 
    sum(a .* f.(xp))
  • przetestuj dla wielomianów coraz większych stopni
  • sprawdz kiedy przestaje być dokładna,
  • podaj związek z twierdzeniem o stopniu dokładności kwadratury Gaussa
In [7]:
function integrate(f, k)
    (xp, w) = gauss(Float64, k)
    sum(w .* f.(xp))  
end
Out[7]:
integrate (generic function with 1 method)
In [8]:
eps = 1.0e-9

for qdeg = 1:5     
    exact = true
    pdeg = 0
    p = Poly([rand()])
    
    while exact   
        p *= Poly([rand(), rand()])
        pdeg += 1
        
        res = polyint(p, -1, 1)
        qres = integrate(p, qdeg)
        diff = abs(qres - res)
        
        print("qdeg = ", qdeg, ", pdeg = ", pdeg)
        println(" | res = ", res, ", qres = ", qres, ", diff = ", diff)
        
        if diff > eps
            exact = false
        end
    end
    
    println("qdeg = ", qdeg, " ==> exact up to ", pdeg - 1, "-degree poly\n")
end

qdeg = 1, pdeg = 1 | res = 0.007227574309893965, qres = 0.007227574309893967, diff = 1.734723475976807e-18
qdeg = 1, pdeg = 2 | res = 0.0265483947625514, qres = 0.002213622556468458, diff = 0.024334772206082943
qdeg = 1 ==> exact up to 1-degree poly

qdeg = 2, pdeg = 1 | res = 0.8846063899545867, qres = 0.8846063899545866, diff = 1.1102230246251565e-16
qdeg = 2, pdeg = 2 | res = 0.6938713861098875, qres = 0.6938713861098874, diff = 1.1102230246251565e-16
qdeg = 2, pdeg = 3 | res = 0.46299023509292503, qres = 0.462990235092925, diff = 5.551115123125783e-17
qdeg = 2, pdeg = 4 | res = 0.3829535561935078, qres = 0.3828213654033437, diff = 0.0001321907901640773
qdeg = 2 ==> exact up to 3-degree poly

qdeg = 3, pdeg = 1 | res = 0.0892009113782127, qres = 0.08920091137821269, diff = 1.3877787807814457e-17
qdeg = 3, pdeg = 2 | res = 0.037908205961984805, qres = 0.03790820596198481, diff = 6.938893903907228e-18
qdeg = 3, pdeg = 3 | res = 0.025296053590456842, qres = 0.025296053590456846, diff = 3.469446951953614e-18
qdeg = 3, pdeg = 4 | res = 0.024313494005645087, qres = 0.02431349400564509, diff = 3.469446951953614e-18
qdeg = 3, pdeg = 5 | res = 0.008178127377284789, qres = 0.00817812737728479, diff = 1.734723475976807e-18
qdeg = 3, pdeg = 6 | res = 0.008294667379781398, qres = 0.008287103928431601, diff = 7.563451349796441e-6
qdeg = 3 ==> exact up to 5-degree poly

qdeg = 4, pdeg = 1 | res = 0.22134387729769606, qres = 0.22134387729769606, diff = 0.0
qdeg = 4, pdeg = 2 | res = 0.23269244923702026, qres = 0.23269244923702032, diff = 5.551115123125783e-17
qdeg = 4, pdeg = 3 | res = 0.1254945893926107, qres = 0.12549458939261077, diff = 5.551115123125783e-17
qdeg = 4, pdeg = 4 | res = 0.12147968265377562, qres = 0.12147968265377565, diff = 2.7755575615628914e-17
qdeg = 4, pdeg = 5 | res = 0.08155106040606931, qres = 0.08155106040606934, diff = 2.7755575615628914e-17
qdeg = 4, pdeg = 6 | res = 0.05929599433400565, qres = 0.059295994334005675, diff = 2.7755575615628914e-17
qdeg = 4, pdeg = 7 | res = 0.07815430782669215, qres = 0.0781543078266922, diff = 4.163336342344337e-17
qdeg = 4, pdeg = 8 | res = 0.07977564337278725, qres = 0.07977307266141405, diff = 2.570711373203749e-6
qdeg = 4 ==> exact up to 7-degree poly

qdeg = 5, pdeg = 1 | res = 1.4055063008467599, qres = 1.4055063008467596, diff = 2.220446049250313e-16
qdeg = 5, pdeg = 2 | res = 1.150899395135486, qres = 1.1508993951354858, diff = 2.220446049250313e-16
qdeg = 5, pdeg = 3 | res = 0.6984474788921667, qres = 0.6984474788921667, diff = 0.0
qdeg = 5, pdeg = 4 | res = 0.5155788979868798, qres = 0.5155788979868798, diff = 0.0
qdeg = 5, pdeg = 5 | res = 0.2591639355508969, qres = 0.2591639355508969, diff = 0.0
qdeg = 5, pdeg = 6 | res = 0.18325819262205012, qres = 0.18325819262205018, diff = 5.551115123125783e-17
qdeg = 5, pdeg = 7 | res = 0.05413414390117652, qres = 0.05413414390117653, diff = 1.3877787807814457e-17
qdeg = 5, pdeg = 8 | res = 0.029574609192085183, qres = 0.02957460919208519, diff = 6.938893903907228e-18
qdeg = 5, pdeg = 9 | res = 0.021075173671651685, qres = 0.021075173671651688, diff = 3.469446951953614e-18
qdeg = 5, pdeg = 10 | res = 0.009244358267261785, qres = 0.009244236696527205, diff = 1.2157073457988443e-7
qdeg = 5 ==> exact up to 9-degree poly

Zadanie 3

Skorzystaj z rozwiązania zadania 2 do napisania funkcji liczącej całki w dowolnym przedziale $\int_{a}^{b} f(x) dx$

dokonując normalizacji do $\int_{-1}^{1} F(z) dz$

podstawiając:

$x=\frac{b+a}{2}+ \frac{b-a}{2} z $ oraz

$dx =\frac{b-a}{2} dz $

Przetestuj działanie na kilku przykładach i sprawdź z wynikami otrzymanymi analitycznie.

In [9]:
function integrate(f, k, a, b)
    fn = x -> f((b + a)/2 + (b - a)/2 * x)
    return (b - a)/2 * integrate(fn, k)
end
Out[9]:
integrate (generic function with 2 methods)
In [10]:
p1 = Poly([1., 2., 3.]) # 3x^2 + 2x + 1 => x^3 + x^2 + x
println("polyint:   ", polyint(p1, -0.4, 3))
println("integrate: ", integrate(p1, 4, -0.4, 3))

polyint:   39.304
integrate: 39.30400000000001
In [11]:
println("polyint:   ", exp(3.3) - exp(-1.1))
println("integrate: ", integrate(exp, 5, -1.1, 3.3))

polyint:   26.779767836959802
integrate: 26.77975216952316

Zadanie 4

Głowną funkcją pakietu QuadGK jest adaptacyjna funkcja guadgk używająca całkowania Gauss-Kronroda

  • użyj tej funkcji do policzenia całki dla przykładowego wielomianu.
  • funkcja ta ma możliwość liczenia również całek do nieskończoności Policz całkę od minus do plus nieskonczonosci standardowego rozkładu normalnego Gaussa $ \frac{1}{\sqrt{2\pi}}exp(\frac{-x^2}{2})$
In [12]:
p1 = Poly([1., 2., 3.]) # 3x^2 + 2x + 1 => x^3 + x^2 + x
println("polyint: ", polyint(p1, -5, 7))
println("quadgk: ", quadgk(p1, -5, 7))

polyint: 504.0
quadgk: (504.0, 1.1368683772161603e-13)
In [13]:
gaussian = x -> (1 / (sqrt(2*pi))) * exp(-(x*x)/2)
println("quadqk for gaussian: ", quadgk(gaussian, -Inf, Inf))
plot(gaussian, -4, 4, label = "gaussian distribution", size = (900, 350))

quadqk for gaussian: (1.0000000000032583, 1.4395584941504537e-8)
Out[13]:
-4 -2 0 2 4 0.0 0.1 0.2 0.3 0.4 gaussian distribution

Zadanie 5

Napisz własną funkcję całkującą metodą prostokątów albo trapezów. Narysuj wykres funkcji błędu w stosunku do wyniku otrzymanego analitycznie w zaleznosci od ilosci potrzebnych przedziałów dla przykładowego wielomianu.

In [14]:
function trapezoid(f, a, b, n)
    sum = 0
    dx = (b - a) / n
    xs = range(a, stop = b, length = n + 1)
    for i = 1:n
        sum = sum + (f(xs[i]) + f(xs[i+1])) * dx / 2
    end
    sum
end
Out[14]:
trapezoid (generic function with 1 method)
In [15]:
p1 = (1/20) * Poly([64., 10., -27., -11., 3., 1.])
println(p1)
a = -4
b = 2
plot(x -> p1(x), a-0.1, b+0.1, label = "p1 - sample poly")

for n = 1:2:7
    xs = range(a, stop = b, length = n + 1)
    p = plot!(xs, p1(xs), label = "trapezoid n = $n")
end
plot!(size = (900, 350), legend = :bottomleft)

Poly(3.2 + 0.5*x - 1.35*x^2 - 0.55*x^3 + 0.15000000000000002*x^4 + 0.05*x^5)
Out[15]:
-4 -3 -2 -1 0 1 2 -1 0 1 2 3 4 5 p1 - sample poly trapezoid n = 1 trapezoid n = 3 trapezoid n = 5 trapezoid n = 7
In [16]:
exact = polyint(p1, a, b)
ns = 1:2:200
err = [abs(exact - trapezoid(p1, a, b, n)) for n in ns]
scatter(ns, err, xlabel = "subintervals", ylabel = "absolute error",
    label = "abs_error(subintervals)", yscale = :log10, size = (900, 400))
Out[16]:
0 50 100 150 200 10 - 3 10 - 2 10 - 1 10 0 10 1 subintervals absolute error abs_error(subintervals)